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ABSTRACT 

We construct a simple dynamical model of the massive globular cluster to Cen. The 
model includes simple treatments of dynamical evolution in a galactic tide, and evolu- 
tion of single stars. Binary stars and rotation are neglected. The model approximately 
fits observational data on the surface brightness profile, the profile of radial velocity 
dispersion, and the main sequence mass function at two radii. 

Key words: stellar dynamics - methods: miscellaneous - stars: luminosity function, 
mass function - globular clusters: individual: uj Cen 



1 INTRODUCTION 

Dynamical modelling of individual globular clusters has a 
long history (see, for example, the review by Meylan & Heg- 
: gie 1997). For the most part the models used have been 
variants of the King-Michie model. While incorporating es- 
sential aspects of stellar dynamics they have nothing to say 
about dynamical evolution, which has been the focus of the- 
ory for many years. 

Dynamical evolutionary models, based on Fokker- 
Planck codes, have been constructed for a number of in- 
dividual clusters by members of Cohn's group (Grabhorn et 
al 1992, Dull et al 1997, Drukier 1995). These models were 
tested against available observational data on the surface 
brightness profiles, individual stellar radial velocities, some- 
times ground-based mass functions, and even data on the 
exotic stellar components in a cluster (Phinney 1993). 

In the meantime a wealth of high quality data on the 
main sequence mass function down to near the H-burning 
limit has been acquired, thanks mainly to the high resolu- 
tion of HST (e.g. Elson et al 1995, Cool et al 1996, Santiago 
et al 1996, von Hippel et al 1996, Piotto et al 1997, King et 
al 1998, Marconi et al 1998, Pulone et al 1999, Paresce & De 
Marchi 2000, De Marchi et al 2000, Andreuzzi et al 2001, 
where we have restricted the credits to one citation per lead 
author). Though static (King- like) models have been con- 
structed which incorporates some of this data (Anderson 
1997, Sosin & King 1997, Saviane et al 1998, Piotto & Zoc- 
cali 1999, in addition to some of the foregoing references), 
we are not aware of any evolutionary models which do so. 

Our aim in this paper is a first step in this direction, 
i.e. to construct a dynamical evolutionary model of a well 
observed, old, galactic globular cluster, so as to represent 
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its surface brightness profile, profile of radial velocity dis- 
persion, and mass function. The object we have chosen for 
this project is uj Cen (NGC 5139). Being massive, it has 
a long relaxation time, and is therefore dynamically rela- 
tively unevolved. Despite the uncertain effects of evolution 
of stars above the turnoff, the mildness of any dynamical 
evolution makes it relatively straightforward to guess ap- 
propriate initial conditions. Nevertheless, some iteration is 
necessary, and for this purpose a fast method of computing 
dynamical evolution is preferred. We have adopted a slightly 
modified version of a Monte Carlo code. 

In the following section we first summarise the obser- 
vational data we have used. Next we outline the manner in 
which the code was adapted for this specific application, and 
how the output was converted for comparison with observa- 
tional data. We then present our model, and discuss briefly 
how the initial conditions were arrived at. The final section 
summarises our conclusions, and discusses the simplifying 
assumptions on which our work is based. 



2 A MODEL OF uj Cen 

2.1 Observational constraints 

uj Cen is a famous object, and has even had its own confer- 
ence recently. In the proceedings (van Leeuwen et al 2002) 
our present understanding of the many sides of uj Cen is 
summarised in detail. For our purposes, however, the obser- 
vational data will be restricted to three sources: (i) a surface 
brightness profile; (ii) the dispersion of radial velocities; and 
(iii) the mass function at two radii. 

For the surface brightness profile we have adopted the 
data in the compilation by Meylan (1987). The surface 
brightness profile consists of both photometric data and val- 
ues derived from star counts, the latter having been nor- 
malised to fit the former in the overlap region. In modelling 
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it might well be better to treat the two kinds of data sepa- 
rately, but we have followed Meylan in treating the surface 
brightness profile as a single data set. The surface brightness 
has been adjusted approximately for interstellar absorption. 
Much doubt was cast on the reliability of this surface bright- 
ness profile by van Leeuwen et al (2000), on the basis of a 
proper motion study, but more recently the results have been 
reconciled (van Leeuwen & Le Poole 2002). 

Our radial velocity data is taken from Meylan et al 
(1995). We have used the binned data presented in this pa- 
per, rather than individual values in Meylan's catalogue. 
This implies that the largest radius is at about 28pc, com- 
pared with a nominal tidal radius of about 70pc (Meylan 
1987). Therefore at large radii the surface brightness profile 
is the only constraint. 

For the mass function we have used the luminosity func- 
tions presented by De Marchi (1999), converted to a mass 
function using the model of Baraffe et al (1997) for metal- 
licity [M/H] — —1.5. We have assumed that the data of De 
Marchi give the numbers of observed stars per magnitude 
bin in the two fields (with his stated adjustment by 0.1 dex 
for the outer field). This assumption is consistent with the 
stated numbers of objects in each field. The two fields for 
which mass functions are given are stated to be at radii of 
about 7' and 4.6', and these are the values we have adopted. 
Nevertheless, from an inspection of the image fields in re- 
lation to lu Cen, we consider that the outer field is centred 
more nearly at 7.4'. We have also corrected these radii for 
the ellipticity, using values from Geyer et al (1983). A more 
recent determination of the ellipticity (van Leeuwen & Le 
Poole 2002) suggests that the ellipticity is even larger. 

2.2 Features of the Monte Carlo code 

The Monte Carlo code (Giersz 1998, 2001) models a spher- 
ically symmetric stellar system by a number of spherical 
shells characterised by radius, energy and angular momen- 
tum. As described in these papers the number of such shells 
equals the number of stars in the real system. In the present 
application, however, the number of shells is much smaller 
than the number of stars, just as in the earlier formulations 
by Henon (1971) and Stodolkiewicz (1982), and so each shell 
is referred to as a superstar. Each superstar represents a 
large number of stars of the same stellar mass, energy and 
angular momentum. In this investigation the initial num- 
ber of shells was chosen in the range from 1024 to 16384. 
Where individual models are discussed below, 16384 shells 
were used. 

Very briefly, the code chooses the radius of each star 
in a manner appropriate to its energy and angular momen- 
tum. The energy and angular momentum are adjusted in 
a manner dictated by the theory of two-body relaxation. 
These processes are repeated until the required evolution 
time (which we assume to be 12Gyr) has elapsed. 

For dynamical purposes, stellar evolution is treated in 
a very simple manner. At the start, each star is assigned an 
evolution time (depending on its initial mass) according to 
the prescription adopted by Chernoff & Weinberg (1990). At 
this time the star is replaced by a degenerate remnant whose 
mass is also determined as in Chernoff & Weinberg, except 
for one point: for an initial mass in the range 4.7Mq < 
m < 8.OM0 the remnant is assumed to be a neutron star 



of 1.4Mq. Neutron stars are given no kick, and so all are 
retained, except those which may escape by relaxation or 
tidal overflow; see also Sec. 3. 3. 

Now we discuss briefly the data output from the Monte- 
Carlo code. At any time a model is defined by the radius, 
energy, angular momentum and mass of its supershells. The 
resulting surface brightness and radial velocity dispersion 
profiles may be very noisy, for two reasons. First, each shell 
gives a cusp in surface brightness at its edge. Secondly, these 
profiles are dominated by the relatively small number of 
evolving stars. It is better, therefore, to represent each su- 
perstar by a space density, corresponding to its orbital mo- 
tion. As a compromise, we represented each superstar by 
100 radii, chosen with the correct radial distribution for a 
superstar with given energy and angular momentum. 

2.3 Conversion to observational data 

Each shell of nominal radius r is taken to represent a uniform 
spherical shell of radii 0.9r and l.lr, to further reduce the 
effects of the cusp which would exist at the projected edge 
of a thin shell. This shell also corresponds to known values 
of the radial and transverse velocity. It is therefore easy to 
compute a density-weighted velocity dispersion along a line 
of sight. Because the observed radial velocities are obtained 
for giants, only shells corresponding to non-degenerate stars 
above O.7M0 were included. (Though not all such stars 
would be giants, it was assumed that mass segregation in 
the range of masses from O.7M0 to the turnoff mass could 
be ignored.) 

The mass functions were obtained in a similar way. 

The most problematic area is creation of the surface 
brightness profile, which requires computation of the lumi- 
nosity of each star from its mass and age (12 Gyr). For 
main sequence stars we used the formulae of Eggleton et al 
(1989), but scaled the evolution time (i.e. the time for the 
end of non-degenerate phases of evolution) to coincide with 
those used in the Monte Carlo code (i.e. those of Chernoff 
& Weinberg 1990). 

Applied to evolving stars, this approach led to the oc- 
currence of a very few shells with very high luminosity, which 
produced a very rough ("bumpy") surface brightness pro- 
file. Therefore we computed the time-averaged luminosity 
during these phases of evolution, and assigned this lumi- 
nosity to each evolving star (i.e. post-main sequence but 
non-degenerate stars). For this purpose we used the code of 
Hurley et al (2000). We checked that the mean luminosity 
of stars brighter than my = 16 in the catalogue of Lynga 
(1996) is approximately consistent with what would be ob- 
tained from the code. 

The surface brightness was corrected for extinction. A 
simple bolometric correction was also applied (Reed 1998)R 



3 FINDING A MODEL OF u> Cen 
3.1 Initial conditions 

Our initial model is a King model (King 1966), specified 
by its total mass, M, and scaled central potential, Wq. The 

1 Note that, in his formula (5), T should be replaced by 10 -4 T 
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initial tidal radius was set from observational estimates of 
current values (M = 3.9 x 1O 6 M [Pryor & Meylan 1993], 
rt = 63.9pc, from data in Trager et al (1993) and Peterson 
(1993)), assuming that r t oc M 1,z . 

Guided by recent observational data (Kroupa 2001) we 
adopted an initial mass function in the form of a continuous 
broken power law 

f(rr,\ nr i m ~ ai ^ H mi < m < m b ; 
11 > \m- a \ if m ft <m<m 2 . 

We quickly realised that ai was quite tightly constrained 
near ct\ = 1 by the mass functions, and adopted this value. 
We also fixed mi = O.1M and m.2 = 15M®. The value of 
mi is a little lower than the lowest mass included in the 
mass functions. Specification of the upper mass limit m.2 is 
relatively unimportant in our models: because the number of 
superstars is so modest, all shells have masses considerably 
below m,2, unless Q2 is rather low. The mass function is 
therefore specified by the mass at the break point between 
the two power laws, m b , and the slope of the initial mass 
function for higher masses, a.i. 

In summary, each initial model is specified by the four 
parameters M, Wo, 0.1, m&. 

3.2 Exploration of initial conditions 

Using 4096 shells, the computation of a single model takes a 
few minutes on a 400MHz Sun workstation. After the Monte 
Carlo code has run, the output can then be compared with 
the three kinds of data, i.e. the profiles of surface bright- 
ness and velocity dispersion, and the mass functions at two 
observed radii. This was done both visually and by a calcu- 
lation of \ 2 - F° r the latter purpose, estimates of the errors of 
the observational data were adopted. Because of the Monte 
Carlo nature of the code, the predictions of each model are 
also subject to statistical uncertainty, but no attempt was 
made to quantify this for purposes of computation of \ 2 ■ 

After preliminary examination of a number of models, 
the parameter space was explored somewhat more system- 
atically, but still manually, by considering the effect of vari- 
ation in each of the parameters. The conclusions are sum- 
marised in Table 1, which also gives the ranges of values of 
the four parameters outside which the fit was observed to 
deteriorate grossly (as judged both by graphical display and 
by the values of \ 2 )- At this stage the best parameter values 
found were approximately M — 1O 7 M , Wo = 8, a.i = 1.9 
and m b — O.6M . 

In order to improve these preliminary values several 
methods were tried, and we describe here the two most suc- 
cessful ones. The first was simply to conduct a Monte Carlo 
search of the range constrained by the values in Table 1, 
i.e. by uniform random sampling of the corresponding hy- 
percube in parameter space. By plotting the resulting val- 
ues of x 2 against each parameter, it was quite easy to de- 
termine the best values with acceptable accuracy. This is 
a relatively slow method, however, and requires thousands 
of Monte Carlo runs. A faster and automatic method, re- 
quiring only of order 50 runs, treats our problem as one of 
stochastic optimization, the stochasticity arising from the 
nature of the Monte Carlo method used for the dynamical 
evolution. Known simply as DIRECT, for "Dividing RECT- 
angles" (Jones 2001), it proceeds by subdividing the search 



domain in a manner that balances global and local searches 
for a minimum. 

3.3 Features of a typical model 

Both methods described in the previous subsection led to 
fairly consistent conclusions. The first method yielded initial 
conditions M ~ 1.0 x 1O 7 M , Wo ~ 7.7, a 2 ~ 1.9 and m b ~ 
O.6M , while DIRECT gave results in the ranges 0.94 x 
1O 7 M < M < 1.23 x 1O 7 M , 7.4 < W < 7.9, 1.95 < 
a 2 < 2.12 and O.63M < m b < 1.14M©. One of the best 
models is illustrated in Fig.l. The current mass is 3.6 x 
1O 6 M , which is perhaps a little too small: the resulting 
tidal radius is perhaps a little too small to account for the 
surface brightness profile at the largest radii. On the other 
hand the modelling of the tide as a cutoff is inaccurate near 
the tidal radius, and so a good fit here may not be achievable. 

We have not tried to adjust the initial mass function 
on the lower main sequence to improve the detailed fit with 
the mass function. Of greater concern is the fact that the 
model mass functions are often slightly but systematically 
too low or high, at the inner and outer radii, respectively. 
The suggestion that the ellipticity exceeds the value we used 
(see Sec. 2.1) would help. 

In attempting to construct a multi-mass King model 
for u) Cen, Meylan (1987) drew attention to the need for 
heavy remnants, by which we mean here both neutron stars 
and white dwarfs. Our models include such remnants, which 
arise from the evolution of stars above the turnoff mass in 
our mass function. Their proportion by mass at the present 
day, in our best models, is of order 50%, and even so it is 
not possible to quite reach the observed velocity dispersion 
at small radii. This problem worsens if all neutron stars are 
removed at birth (i.e. it is assumed that each receives a 
kick exceeding the escape speed). Though the total fraction 
of degenerate stars declines only to about 40%, the central 
velocity dispersion drops to about 12km/s. 

Several other features of these models may be of in- 
terest. The models are mildly anisotropic. If the anisotropy 
parameter (3 is defined by (3 = 1 — (v% )/(«?), where v r , v t are 
the radial and transverse velocities in the plane of the sky 
(i.e. as measured by proper motions), we find that (3 varies 
from a value of about within the innermost 2pc to about 
—0.15 at a radius of 10 pc, and then rises towards as the 
tidal boundary is approached. At 20 pc, close to the radius 
where King & Anderson (2002) found only mild anisotropy, 
13 ~ -0.05. 

The models also exhibit mild mass segregation. In the 
model exhibited in Fig.l, the mean mass of unevolved stars 
is nearly O.4OM at all (projected) radii less than about 
lpc, and about 0.34 beyond lOpc; the mean mass declines 
steadily between 1 and lOpc. The primordial value (over 
the same range of stellar masses) was O.35M . For a single 
power law f(m)dm oc m~ a dm between the minimum mass 
and turnoff, the variation of mean mass corresponds to a 
variation of a of about 0.5. Relative to the centre, therefore, 
there is an excess of stars of lowest mass at the outside of 
the cluster by a factor of order 3 (0.5dex). Anderson (2002) 
has observed slight mass segregation by comparing the lumi- 
nosity function at the centre relative to a field at about 7'. 
His faintest stars are underabundant at the centre by about 
0.2dex, but do not extend to such low masses. 
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Table 1. Preliminary parameter values. 



Parameter and range 



effect of increase 



effect of decrease 



6 x 10 6 M Q < M < 1.4 x 10 7 
7.1 < Wo < 8.6 
1.7 <a 2 < 2.2 



O.6M < m b < 0.95M Q 



mf and v 2 too great 
sfb too concentrated 
mf and sfb too high 



sfb too concentrated, 



mf and v 2 too low 
underluminous at centre 
model underluminous; 
central v 2 too small; 
mf at inner radius poor 
centre overluminous, mf wrong 
mf at smaller radius too low 



Notes: mf = mass function; sfb = surface brightness profile; v = radial velocity dispersion 
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Figure 1. One of the best models. Here M = 0.94 X 1O 7 M0, Wo = 7.6, 02 = 1-95 and rn b = 
O.63M0; the number of superstars was 16384. Left: logarithmic surface brightness profile (in 
units of 10.00 V mag per square arc minute) against R (pc); lower right: velocity dispersion 
profile (km/s); upper right: mass function (stars per unit mass per square arc minute; the upper 
and lower plots correspond to the inner and outer observed radii, respectively). 



The fact that the signatures of dynamical evolution in 
lo Cen are not great is no surprise, but this is perhaps the 
first time they have been quantified theoretically. 



4 DISCUSSION AND CONCLUSIONS 
4.1 Discussion 

Before summarising the tentative findings of this study, it is 
important to set out some aspects of u> Cen which we have 
not included. 

The most important constituent we have omitted is any 
population of binary stars. At the expense of introducing 
further parameters, it would have been possible to do so 
by treating the population as dynamically "inert", i.e. as 
simply a population of slightly more massive stellar objects. 
Nevertheless it would have been desirable to treat their stel- 
lar evolution in a fundamentally different way from that of 
the single stars, and so for this exploratory study they were 
neglected entirely. If included, they might have helped to 
deepen the potential well and increase the central velocity 
dispersion, without unduly distorting the surface brightness 
there. They might also assist the retention of neutron stars, 



with the same result. Our treatment of the evolution of sin- 
gle stars could also be improved significantly. 

Dynamically, we have taken no notice of the fact that 
ui Cen is rotating (see Merritt et al 1997). For purposes of 
dynamical evolution it seems that rotation does not play a 
dominant role (Spurzem 2001), but our practical reason for 
omitting rotation is that the code cannot cope with it. The 
rotation of ui Cen may be a symptom of past mergers (Norris 
et al 1997), and this possibility is ignored here also. 

In fitting models to observational data we have ignored 
the kinematical evidence from internal proper motions (van 
Leeuwen et al 2000). This data appears to be entirely consis- 
tent with the radial velocity, which we did employ, and also 
covers a similar range of radius within the cluster. Another 
reason for neglecting this data is that we have it in mind to 
apply our methods to several other clusters for which such 
data is lacking entirely. 

A significant dynamical mechanism that we have also 
ignored is the time taken for stars to escape. As Baumgardt 
(2001) has shown, the effect is that the lifetime in a tidal 
field is not proportional to the relaxation time, as we would 
find using our Monte Carlo code. On the other hand we have 
also simplified the treatment of the tidal boundary condition 
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by supposing that the tide is steady, as for a cluster on a 
circular galactic orbit. We have made no attempt to model 
the tidal debris of u) Cen, whose mass is considerable (Leon 
et al 2000). 

An issue which we would have liked to address is the 
uniqueness of the model initial conditions that we have 
found. While we have explored the parameter space in var- 
ious ways, and have tried to place bounds on the parame- 
ters which give acceptable models, perhaps radically differ- 
ent models are possible. 

Our best model is no more than a tolerable fit to the 
data with which it has been compared. On the other hand 
the data itself is not without problems, such as the diffi- 
culty of converting from a magnitude distribution to a mass 
function. 



4.2 Conclusions 

We have found that the surface brightness profile, velocity 
dispersion profile and mass function of uo Cen can be fitted 
approximately by the dynamical evolution, over 12Gyr, of a 
cluster with the following initial conditions: the initial mass 
is about 1.1 x 1O 7 M0, the initial tidal radius is about 90pc, 
and the initial model is a King model with a scaled central 
potential Wo = 7.7 approximately; the initial mass function 
is a broken power law, with slopes of about 1 and 1.9 re- 
spectively below and above the break-point mass of about 
O.6M . 

The resulting present mass and tidal radius are about 
3.6 x lO 6 A'/0 and 61pc, respectively. The current proportion 
of mass in heavy remnants in our model is about 55%. 
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